ESM 244 Lab 1

# Reading in our trees data (use read_csv for more options - better than read.csv)
sf_trees <- read_csv(here('data', 'sf_trees', "sf_trees.csv"), 
                     show_col_types = FALSE)

# Always inspect the data first 
summary(sf_trees)

##     tree_id       legal_status         species            address         
##  Min.   :     1   Length:192987      Length:192987      Length:192987     
##  1st Qu.: 52602   Class :character   Class :character   Class :character  
##  Median :120862   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :126529                                                           
##  3rd Qu.:202608                                                           
##  Max.   :261546                                                           
##                                                                           
##    site_order       site_info          caretaker              date           
##  Min.   :-50.000   Length:192987      Length:192987      Min.   :1955-09-19  
##  1st Qu.:  1.000   Class :character   Class :character   1st Qu.:1995-01-09  
##  Median :  2.000   Mode  :character   Mode  :character   Median :2001-06-12  
##  Mean   :  4.579                                         Mean   :2000-10-12  
##  3rd Qu.:  4.000                                         3rd Qu.:2008-09-25  
##  Max.   :501.000                                         Max.   :2020-01-25  
##  NA's   :1634                                            NA's   :124610      
##       dbh           plot_size            latitude       longitude     
##  Min.   :   0.00   Length:192987      Min.   :37.51   Min.   :-138.3  
##  1st Qu.:   3.00   Class :character   1st Qu.:37.74   1st Qu.:-122.5  
##  Median :   7.00   Mode  :character   Median :37.76   Median :-122.4  
##  Mean   :   9.95                      Mean   :37.77   Mean   :-122.4  
##  3rd Qu.:  12.00                      3rd Qu.:37.78   3rd Qu.:-122.4  
##  Max.   :9999.00                      Max.   :47.27   Max.   :-122.4  
##  NA's   :41819                        NA's   :2832    NA's   :2832

names(sf_trees)

##  [1] "tree_id"      "legal_status" "species"      "address"      "site_order"  
##  [6] "site_info"    "caretaker"    "date"         "dbh"          "plot_size"   
## [11] "latitude"     "longitude"

# Viewing the top 5 most common legal statuses and their respective tree count numbers 
top_5_status <- sf_trees %>% 
  group_by(legal_status) %>% 
  summarize(tree_count = n()) %>% 
  slice_max(tree_count, n = 5) %>% # just keeps the maximum/largest values of the number of top observations we want 
  arrange(-tree_count) # arranges selected choices from greatest to least 

Creating a GG Plot of Top 5 SF trees legal status data

ggplot(data = top_5_status, aes(x = fct_reorder(legal_status, tree_count), y = tree_count)) + 
       # fct_reorder - reorders legal statuses from default alphabetical order to lowest - the greatest tree count
  geom_col() + 
  labs(y= 'Tree Count', x = 'Legal Status') + 
  coord_flip() + # flips the axes 
  theme_minimal()

Keep obs where legal status is “Permitted Site” and caretaker is “MTA”

# sf_trees$legal_status %>% unique()
permitted_mta <- sf_trees %>% 
  filter(legal_status == "Permitted Site" & caretaker =="MTA")

Keep Blackwood acacia trees, only legal_status, date, latitude, longitude columns

blackwood_acacia <- sf_trees %>% 
  filter(str_detect(species, "Blackwood Acacia")) %>% 
  select(legal_status, date, latitude, longitude)

# Make a plot of SF acacia 
ggplot(data = blackwood_acacia, aes(x = longitude, y = latitude)) +
  geom_point()

Use the tidyr::separate() to separate one column into multiple columns, and tidyr::unite() to rejoin

# Using separate to break column into species and scientific name 

sf_trees_sep <- sf_trees %>% 
  separate(species, into = c('spp_scientific', 'spp_common'), sep = ' :: ')

sf_trees_unite <- sf_trees %>% 
  unite('id_status', tree_id:species) # or could use commas to join specific columns not in a row 

Make some actual maps!

# Spacial / simple features map 
blackwood_acacia_sf <- blackwood_acacia %>% 
  drop_na(longitude, latitude) %>% 
  st_as_sf(coords = c('longitude', 'latitude')) # from the sf package - turn spacial & temporal data and turn it into a geometry feature (a point) 

st_crs(blackwood_acacia_sf) <- 4326
# telling R this coordinate reference system is important
# numeric code - indicates basic latitude and longitude coordinate reference system in degrees 

ggplot(data = blackwood_acacia_sf) + 
  geom_sf(color = 'darkgreen') + # plots geometrical special features (always looks for the geometry column as the aesthetics to map the x & y)
  theme_minimal()

Read in SF streets data

sf_map_sf <- read_sf(here('data', 'sf_map','tl_2017_06075_roads.shp')) %>%  # reads in spacial data as a simple feature 

# st_crs(sf_map_sf) # tells us about the crs 
# want to change the crs (before we knew what it was so we could just type it in, but this one is different)

# we want to tranform it into the one we want
st_transform(4326)

st_crs(sf_map_sf) # now it's where we want it to be 

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

Let’s create a plot now

ggplot() +
  geom_sf(data = sf_map_sf, size = 0.1, color = 'darkgrey') + 
  geom_sf(data = blackwood_acacia_sf, color = 
          'red', size = 0.5) + 
            theme_void() + 
            labs(title = "Blackwood acacias in San Francisco")

Interactive map!

tmap_mode('view')
tm_shape(blackwood_acacia_sf) +
  tm_dots()

# creating a shape using this data with dot shapes